import java.io.File;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.Collections;
import java.util.InputMismatchException;
import java.util.Iterator;
import java.util.LinkedHashMap;
import java.util.LinkedList;
import java.util.ListIterator;
import java.util.Map;
import java.util.Scanner;
import java.util.TreeMap;

import javax.xml.parsers.DocumentBuilderFactory;
import javax.xml.parsers.DocumentBuilder;

import org.w3c.dom.Document;
import org.w3c.dom.Element;
import org.w3c.dom.NodeList;
import org.w3c.dom.Node;

import org.biojava3.core.sequence.AccessionID;
import org.biojava3.core.sequence.DNASequence;
import org.biojava3.core.sequence.DataSource;
import org.biojava3.core.sequence.RNASequence;
import org.biojava3.core.sequence.TaxonomyID;
import org.biojava3.core.sequence.compound.AmbiguityDNACompoundSet;
import org.biojava3.core.sequence.compound.NucleotideCompound;
import org.biojava3.core.sequence.storage.ArrayListSequenceReader;
import org.biojava3.core.sequence.compound.AmbiguityRNACompoundSet;
import org.forester.io.parsers.phyloxml.PhyloXmlParser;
import org.forester.phylogeny.Phylogeny;
import org.forester.phylogeny.PhylogenyMethods;
import org.forester.phylogeny.PhylogenyNode;
import org.forester.phylogeny.iterators.PhylogenyNodeIterator;

public class OpenFasta {
	private LinkedList<Site> allSequenceSites = new LinkedList<Site>();
	private LinkedList<RNASequence> xmlSequences = new LinkedList<RNASequence>();
	private LinkedList<RNASequence> nonRedundantSequences = new LinkedList<RNASequence>();
	private LinkedList<RNASequence> ampliconSequencesRNA = new LinkedList<RNASequence>();
	private LinkedList<DNASequence> ampliconSequences = new LinkedList<DNASequence>();
	private LinkedList<String> rankList = new LinkedList<String>();
	private LinkedList<String> silvaRankList = new LinkedList<String>();
	private LinkedHashMap<String, Integer> classifiedRanks = new LinkedHashMap<String, Integer>();
	private LinkedHashMap<Integer, LinkedList<String>> lowerRankList = new LinkedHashMap<Integer, LinkedList<String>>();
	private LinkedList<Phylum> chosenRankList = new LinkedList<Phylum>();
	private LinkedList<RNASequence> phylumAmplicons = new LinkedList<RNASequence>();
	private LinkedHashMap<String, String> subphylaList =  new LinkedHashMap<String, String>();
	private LinkedList<String> databaseOriginale = new LinkedList<String>();
	private LinkedList<String> databaseReformatio = new LinkedList<String>();
	private Phylogeny[] phylogenies = null;
	private NCBITaxonomy[] ncbiTaxonomies = null;
	private String currentRank = null; 
	private String[] samples = {"Sample01", "Sample02", "Sample03", "Sample04", "Sample05", "Sample06", "Sample07", "Sample08", "Sample09", "Sample10"};

	private String referenceSilvaName = "ShyCer61";
	private RNASequence referenceSequence = null;
	private String forwardPrimer = "CAGCAGCCGCGGUAAUUCC";
	private String reversePrimer = "CCCGUGUUGAGUCAAAUUAAGC";
	private BuildPhyloXMLTree mainPhyloTree = null;
	private BuildPhyloXMLTree nonRedundantPhyloTree = null;
	private BuildPhyloXMLTree notFoundTree = null;
//	String referenceSilvaName = "EschC385";
//	static String forwardPrimer = "CAGCAGCCGCGGUAAUAC";
//	static String reversePrimer = "UCGAAUUAAACCACAU";
	private String searchPhylum = null;
	private String xmlOutputFile = "";
	private String ampliconsFile = null;
	private String phyloXMLFile = "";
	private String fileAmpliconsRNA = "";
	private String fileAmpliconsDNASense = "amplicons_DNA_w_" + forwardPrimer + "_" + reversePrimer;
	private String fileAmpliconsDNANoSense = "amplicons_DNA_wo_" + forwardPrimer + "_" + reversePrimer;
	private String fileAmpliconsLength = "";
	private String fileAmpliconsNames = "";
	private String fileSilva2NCBI = "silva2ncbi.map";
	private String rankFile = null;
	private String xmlFile = "";
	
	private String[] rankFiles = {"z_ncbi_kingdom", "z_ncbi_phylum", "z_ncbi_class", "z_ncbi_order", "z_ncbi_family", "z_ncbi_genus"};
	
	public void initialize(String[] arguments) {
		long start = 0;
		long end = 0;
		
		System.out.println("Hello! Welcome to AmpliconGenerator 0.1!");
		if (arguments.length > 0) {
			System.out.println("Now I'm going to read the provided xml file.");
			xmlFile = arguments[0];
			start = System.currentTimeMillis();
			readXMLFile(xmlFile);
			end = System.currentTimeMillis();
			System.out.println("XML file is read! There are " + xmlSequences.size() + " sequences. " + (double)(end-start)/1000 + " s");
		}
		
		int choice = 0;
		while (choice != -1) {
			System.out.println();
			System.out.println("1. Build datastructure of the sequence alignment columns.");
			System.out.println("2. Find the reference sequence in the datastructure.");
			System.out.println("3. Generate amplicons with current primers.");
			System.out.println("4. Choose new primers.");
			System.out.println("5. Print a PhyloXML tree of all sequences to a file.");
			System.out.println("6. Print all data.");
			System.out.println("7. Print silva to NCBI name conversion map using the current sequences.");
			System.out.println("8. Place all non-redundant sequences into a data structure.");
			System.out.println("9. Check the current primers on the current reference sequence.");
			System.out.println("10. Print RNA amplicons to a file.");
			System.out.println("11. Print DNA amplicons to a file.");
			System.out.println("12. Print a PhyloXML tree of nonredundant sequences to a file.");
			System.out.println("13. Build a forester Phylogeny structure from a PhyloXML file.");
			System.out.println("14. Read phyla files.");
			System.out.println("15. Print statistics for phyla.");
			System.out.println("16. Create subtrees of the phyla and print corresponding amplicons to files.");
			System.out.println("17. Create subtree of a phylum and print the statistics.");
			System.out.println("18. Write the taxonomies of the redundant sequences to a file.");
			System.out.println("19. Read NCBI taxonomy files.");
			System.out.println("20. Find rank occurrencies in silva.");
			System.out.println("21. Find all rank occurrences in NCBI.");
			System.out.println("22. Find all rank occurrencies in NCBI for a superkingdom.");
			System.out.println("23. Read taxonomy from classification results.");
			System.out.println("24. Compare two versions of Silva database using fasta files.");
			System.out.println("25. Build a PhyloXML tree from a list of taxonomies.");
			System.out.println("26. Print all sequences to a fasta file.");
			System.out.println("27. Print names and taxa of nonredundant sequences.");
			System.out.println("28. Subsample sediment reads.");
			System.out.println("29. Check non-redundant sequences.");
			System.out.println("30. Generate sequences in specific region.");
			System.out.println("31. Print statistics for the ranks.");
			System.out.println("32. Compare the taxonomies.");
			
			System.out.print("Please select a choice. To quit press -1: ");
			Scanner input = new Scanner(System.in);
			try {
				choice = input.nextInt();
			} catch (InputMismatchException e) {
				e.printStackTrace();
			}
			switch (choice) {
				case 0:	if (!xmlSequences.isEmpty()) 
							System.out.println("Sequences are allready in the datastructure. Please continue.");
						else {
							System.out.println("Write the name of the xml file containing sequence information: ");
							input.nextLine();
							xmlFile = input.nextLine();
							start = System.currentTimeMillis();
							readXMLFile(xmlFile);
							end = System.currentTimeMillis();
							System.out.println("XML file is read! There are " + xmlSequences.size() + " sequences. " + (double)(end-start)/1000 + " s");
						}
						break;
				case 1:	if (!allSequenceSites.isEmpty())
							System.out.println("Data structure containing all alignment columns allready exists!\nIt contains "
									+ allSequenceSites.size() + " sites.");
						else {
							System.out.println("Generating column data from the sequence alignemnt...");
							start = System.currentTimeMillis();
							constructSequencesSites();
							end = System.currentTimeMillis();
							System.out.println("Column data from all aligned sequences is constructed! " + (double)(end-start)/1000 + " s");
							System.out.println("The sequence alignment have " + allSequenceSites.size() + " columns.");
						}
						break;
				case 2:	System.out.println("Searching for the reference sequence " + referenceSilvaName + "...");
						start = System.currentTimeMillis();
						referenceSequence = searchSequenceBySilvaName();
						end = System.currentTimeMillis();
						if (referenceSequence != null)
							System.out.println("Reference sequence " + referenceSequence.getTaxonomy().getSilvaName() + " is found! " + (double)(end-start)/1000 + " s");
						else {
							System.out.println("There is no sequence with silva name " + referenceSilvaName);
							System.exit(0);
						}
						break;
				case 3:	if (!ampliconSequencesRNA.isEmpty()) {
							RNASequence firstRNA = ampliconSequencesRNA.getFirst();
							String beginString = firstRNA.getSequenceAsString().substring(0, 20);
							String endString = firstRNA.getSequenceAsString().substring(firstRNA.getLength()-20);
							System.out.println("Amplicons allready exist and contain sequence starting with " + beginString + " and ending with " + endString);
							System.out.println("Do you still want to generate new amplicons using current primers " + forwardPrimer + " and " + reversePrimer
									+ " ?\nAnswer y/n: ");
							input.nextLine();
							String answer = input.nextLine();
							if (answer.equalsIgnoreCase("n"))
								break;
						}
						System.out.println("Do you want to avoid mismatches? Answer y/n: ");
						input.nextLine();
						String answer = input.nextLine();
						boolean avoidMismatch;
						if (answer.equalsIgnoreCase("y"))
							avoidMismatch = true;
						else
							avoidMismatch = false;
						System.out.println("Generating amplicons using primers " + forwardPrimer + " and " + reversePrimer + "...");
						start = System.currentTimeMillis();
						int numberAmplicons = generateAmplicons(avoidMismatch);
						end = System.currentTimeMillis();
						System.out.println("Found " + numberAmplicons + " amplicons! " + (double)(end-start)/1000 + " s");
						break;
				case 4:	System.out.println("Your current primers are:");
						System.out.println("Forward primer:\t" + forwardPrimer);
						System.out.println("Reverse primer:\t" + reversePrimer);
						System.out.println("Write your new forward primer: ");
						input.nextLine();
						forwardPrimer = input.nextLine();
						System.out.println("Write your new reverse primer: ");
						reversePrimer = input.nextLine();
						System.out.println("Your current primers are:");
						System.out.println("Forward primer:\t" + forwardPrimer);
						System.out.println("Reverse primer:\t" + reversePrimer);
						break;
				case 5:	if (mainPhyloTree != null)
							System.out.println("A XML phylogenetic tree of all organisms is allready printed to " + xmlOutputFile);
						else {
							System.out.println("Printing all organisms into a PhyloXML tree...");
							start = System.currentTimeMillis();
							xmlOutputFile = "allEukaryotesTaxTree.xml";
							mainPhyloTree = new BuildPhyloXMLTree(xmlSequences, xmlOutputFile);
							mainPhyloTree.buildTree();
							end = System.currentTimeMillis();
							System.out.println("PhyloXML tree built in " + (double)(end-start)/1000 + " s");
						}
						break;
				case 6:	start = System.currentTimeMillis();
						if (allSequenceSites.isEmpty())
							System.out.println("Datastructure containing alignment columns is empty! Please choose nr. 1 first.");
						else if (referenceSequence == null)
							System.out.println("Reference sequence does not exist! Please choose nr. 2 first.");
						else {
							System.out.println("Printing all statistics to files...");
							PrintData printAll = new PrintData(allSequenceSites, referenceSequence);
							printAll.printEverything();
							end = System.currentTimeMillis();
							System.out.println("All statistics printed in " + (double)(end-start)/1000 + " s");
						}
						break;
				case 7:	System.out.println("Generating a file containing Silva to NCBI name conversion map...");
						start = System.currentTimeMillis();
						printSilva2NCBIMap(fileSilva2NCBI);
						end = System.currentTimeMillis();
						System.out.println("Silva to NCBI map for MEGAN printed to file!" + (double)(end-start)/1000);
						break;
				case 8:	if (!nonRedundantSequences.isEmpty()) {
							System.out.println("The datastructure allready contains " + nonRedundantSequences.size() + " sequences!\n" +
									"Dou you still want to fill the structure with other amplicons?\n" +
									"Answer y/n: ");
							input.nextLine();
							answer = input.nextLine();
							if (answer.equalsIgnoreCase("n"))
								break;
						}
						else {
							input.nextLine();
						}
						System.out.println("Please choose the file containing the amplicons:");
						ampliconsFile = input.nextLine();
						System.out.println("Building a datastructure containing all non-redundant sequences...");
						start = System.currentTimeMillis();
						buildNonRedundantSequences(ampliconsFile);
						end = System.currentTimeMillis();
						System.out.println("All non redundant sequences placed into a datastructure! " + (double)(end-start)/1000 + " s");
						break;
				case 9:	if (referenceSequence == null)
							System.out.println("Reference sequence does not exist! Please choose nr. 2 first.");
						else {
							System.out.println("Searching for the current primers in the reference sequence.");
							System.out.println("Searching for the forward primer...");
							start = System.currentTimeMillis();
							int indexForward = searchForwardPrimerPosition(forwardPrimer, referenceSequence.getSequenceAsString(), true);
							end = System.currentTimeMillis();
							if (indexForward < 0) {
								System.out.println("Forward primer not found! " + (double)(end-start)/1000 + " s");
								System.exit(0);
							}
							else
								System.out.println("Forward primer found at position " + indexForward + " " + (double)(end-start)/1000 + " s");
							System.out.println(referenceSilvaName + " has " + referenceSequence.getSequenceAsString().length() + " bases.");
							
							System.out.println("Searching for the reverse primer...");
							start = System.currentTimeMillis();
							int indexReverse = searchReversePrimerPosition(reversePrimer, referenceSequence.getSequenceAsString(), true);
							end = System.currentTimeMillis();
							if (indexReverse < 0) {
								System.out.println("Reverse primer not found! " + (double)(end-start)/1000 + " s");
								System.exit(0);
							}
							else
								System.out.println("Reverse primer found at position " + indexReverse + " " + (double)(end-start)/1000 +  " s");
						}
						break;
				case 10:if (ampliconSequencesRNA.isEmpty())
							System.out.println("No RNA amplicons have been generated yet! Please choose nr. 3 first.");
						else {
							fileAmpliconsRNA = "amplicons_RNA_" + forwardPrimer + "_" + reversePrimer;
							fileAmpliconsLength = "amplicons_length_" + forwardPrimer + "_" + reversePrimer;
							fileAmpliconsNames = "amplicons_names_" + forwardPrimer + "_" + reversePrimer;
							System.out.println("Printing RNA amplicon sequences to a file...");
							start = System.currentTimeMillis();
							printAmplicons(fileAmpliconsRNA, SequenceType.RNA, true);
							end = System.currentTimeMillis();
							System.out.println("All different RNA amplicons printed in " + (double)(end-start)/1000 +  " s");
						}
						break;
				case 11:if (ampliconSequences.isEmpty())
							System.out.println("No DNA amplicons have been generated yet! Please choose nr. 3 first.");
						else {
							System.out.println("Printing DNA amplicon sequences to a file...");
							start= System.currentTimeMillis();
							printAmplicons(fileAmpliconsDNASense, SequenceType.DNA, true);
							printAmplicons(fileAmpliconsDNANoSense, SequenceType.DNA, false);
							end = System.currentTimeMillis();
							System.out.println("All different DNA amplicons printed in " + (double)(end-start)/1000 +  " s");
						}
						break;
				case 12:if (nonRedundantSequences.isEmpty())
							System.out.println("The structure contains no sequences! Please choose nr. 8 first.");
						else {
							System.out.println("Printing all organisms with nonredundant sequences into a PhyloXML tree...");
							xmlOutputFile = ampliconsFile + ".xml";
							start = System.currentTimeMillis();
							nonRedundantPhyloTree = new BuildPhyloXMLTree(nonRedundantSequences, xmlOutputFile);
							nonRedundantPhyloTree.buildTree();
							end = System.currentTimeMillis();
							System.out.println("PhyloXML file of organisms from non redundant sequences is generated! " + (double)(end-start)/1000 + " s");
						}
						break;
				case 13:System.out.println("Write the name of the PhyloXML file containing the phylogeny tree: ");
						input.nextLine();
						phyloXMLFile = input.nextLine();
						start = System.currentTimeMillis();
						readPhyloXMLFile(phyloXMLFile);
						end = System.currentTimeMillis();
						System.out.println("The file is parsed into a Phylogeny structure! " + (double)(end-start)/1000 + " s");
						break;
				case 14:rankList.clear();
						input.nextLine();
						rankFile = input.nextLine();
						currentRank = rankFile.substring(rankFile.lastIndexOf("_")+1);
						start = System.currentTimeMillis();
						boolean allOk = readPhylaFiles(rankFile, false);
						end = System.currentTimeMillis();
						if (allOk)
							System.out.println("Ranks are read into the datastructure!" + (double)(end-start)/1000 + " s");
						else
							System.out.println("The files were not read!");				
						break;
				case 15:if (rankList.isEmpty())
							System.out.println("Rank datastructures are empty! Please choose nr. 14 first.");
						else {
							int begin = rankFile.lastIndexOf("_");
							String currentRank = rankFile.substring(begin);
							String phylaFile = phyloXMLFile + "_" + currentRank;
							System.out.println("Printing " + currentRank + " statistics...");
							printPhylaStatistics(phylaFile, false);
						}
						break;
				case 16:if (phylogenies == null) 
							System.out.println("There is no forester Phylogeny structure! Please choose nr. 13 first.");
						else if (chosenRankList.isEmpty()) 
							System.out.println("There exists no list of phyla! Please choose nr. 14 first.");
						else if (ampliconSequencesRNA.isEmpty()) 
							System.out.println("No amplicons have been generated! Please choose nr. 3 first.");
						else {
							start = System.currentTimeMillis();
							createAndPrintAllPhylumAmplicons();
							end = System.currentTimeMillis();
							System.out.println("Amplicons for all phyla created and printed to files! " + (double)(end-start)/1000 + " s");
						}
						break;
				case 17:if (phylogenies == null) 
							System.out.println("There is no forester Phylogeny structure! Please choose nr. 13 first.");
						else if (chosenRankList.isEmpty()) 
							System.out.println("There exists no list of phyla! Please choose nr. 14 first.");
						else if (ampliconSequencesRNA.isEmpty()) 
							System.out.println("No amplicons have been generated! Please choose nr. 3 first.");
						else {
							input.nextLine();
							System.out.println("Write phylum that you want to search: ");
//							String searchPhylum = input.nextLine();
//							createSubtreeAndPrintPhylumStatistics(searchPhylum);
						}
						break;
				case 18: if (nonRedundantSequences.isEmpty()) {
							System.out.println("The structure contains no sequences! Please choose nr. 8 first.");
						}
						findTaxonomy();
						break;
				case 19:start = System.currentTimeMillis();
						readNCBITaxonomy();
						end = System.currentTimeMillis();
						System.out.println("NCBI taxonomy file is read! " + (double)(end-start)/1000 + " s");
						break;
				case 20:if (xmlSequences.isEmpty()) 
							System.out.println("Sequences from Silva are not read! Please choose nr. 0 first");
						else if(ncbiTaxonomies[2] == null) 
							System.out.println("NCBI taxonomy file is not read! Please choose nr. 19 first.");
						else {
							input.nextLine();
							System.out.println("Write rank you want to search for: ");
							String searchRank = input.nextLine();
							start = System.currentTimeMillis();
							findRanksInSilva(searchRank);
							end = System.currentTimeMillis();
							System.out.println("All " + searchRank + " occurrencies found!");
							System.out.println("There are " + silvaRankList.size() + " different " + searchRank + ".");
							System.out.println("There are " + lowerRankList.size() + " additional occurrencies of lower groups.");
						}
						break;
				case 21:if (ncbiTaxonomies[2] == null)
							System.out.println("NCBI taxonomy file is not read! Please choose nr. 19 first.");
						else {
							input.nextLine();
							System.out.println("Write rank you want to search for: ");
							String searchRank = input.nextLine();
							start = System.currentTimeMillis();
							findRanksInNCBI(searchRank, true, 0);
							end = System.currentTimeMillis();
							System.out.println("All " + searchRank + " occurrencies found! " + (double)(end-start)/1000 + " s");
						}
						break;
				case 22:if (ncbiTaxonomies[2] == null) 
							System.out.println("NCBI taxonomy file is not read! Please choose nr. 19 first.");
						else {
							input.nextLine();
							System.out.println("Write rank you want to search for: ");
							String searchRank = input.nextLine();
							System.out.println("Write UID of the superkingdom you want to keep: ");
							String retainUIDString = input.nextLine();
							Integer retainUID = new Integer(retainUIDString);
							start = System.currentTimeMillis();
							findRanksInNCBI(searchRank, false, retainUID);
							end = System.currentTimeMillis();
							System.out.println("All " + searchRank + " occurrencies found! " + (double)(end-start)/1000 + " s");
						}
						break;
				case 23:System.out.println("Current rank is " + currentRank);
						for (String currSample : samples) {
							readClassifiedFasta(currSample + "_.fasta", currentRank);
						}
						break;
				case 24:input.nextLine();
						System.out.println("Write the name of the fasta file containing the original database: ");
						String originalDatabase = input.nextLine();
						System.out.println("Write the name of the fasta file containing the changed database: ");
						String changedDatabase = input.nextLine();
						start = System.currentTimeMillis();
						compareTwoSilvaDatabases(originalDatabase, changedDatabase);
						end = System.currentTimeMillis();
						System.out.println("The two databases are compared and written to files! " + (double)(end-start)/1000 + " s");
						break;
				case 25:input.nextLine();
						System.out.println("Write the name of the file containing taxonomies: ");
						String inFile = input.nextLine();
						System.out.println("Write the name of the file containing the PhyloXML tree: ");
						String outFile = input.nextLine();
						notFoundTree = new BuildPhyloXMLTree(inFile, outFile);
						notFoundTree.buildTreeFromTaxonomy();
						break;
				case 26:start = System.currentTimeMillis();
						printAllToFasta();
						end = System.currentTimeMillis();
						break;
				case 27:start = System.currentTimeMillis();
						printNamesAndTaxa();
						end = System.currentTimeMillis();
						break;
				case 28:input.nextLine();
						System.out.println("Write the name of the file containing flowgrams: ");
	//					String flowFile = input.nextLine();
						System.out.println("Write the name of the fasta file containing sequences: ");
	//					String fastaFile = input.nextLine();
						System.out.println("How many subsamples do you want to generate: ");
	//					String numberSubSamples = input.nextLine();
						start = System.currentTimeMillis();
	//					subsampleSedimentReads(flowFile, fastaFile, numberSubSamples);
						end =System.currentTimeMillis();
						break;
				case 29:if (nonRedundantSequences.isEmpty()) {
							System.out.println("There are no non redundant sequences stored! Please choose nr. 8 first!");
						} else {
							checkNonRedundant();
						}
						break;
				case 30:if(xmlSequences.isEmpty()) {
							System.out.println("There are no sequences in the main database! Please choose nr. 0 first!");
						} else {
							input.nextLine();
							System.out.println("Write the string to associate with the files: ");
							String stringMark = input.nextLine();
							System.out.println("Write the first position: ");
							String firstPosString = input.nextLine();
							Integer startInt = new Integer(firstPosString);
							System.out.println("Write the last position: ");
							String lastPosString = input.nextLine();
							Integer stopInt = new Integer(lastPosString);
							generateAllSequencesInSpecificRegion(stringMark, startInt, stopInt);
						}
						break;
				case 31:printStatisticsForAllRanks();
						break;
				case 32:if(xmlSequences.isEmpty()) {
							System.out.println("There are no sequences in the main database! Please choose nr. 0 first!");
						} else {
							input.nextLine();
							System.out.println("Write the string to search for in the taxonomies: ");
							String taxSearch = input.nextLine();
							compareTaxonomies(taxSearch);
						}
						break;
				default:break;
			}
		}
	}
	
	public void compareTaxonomies(String taxSearch) {
		ListIterator<RNASequence> seqIter = xmlSequences.listIterator();
		int silvaCount = 0;
		int emblCount = 0;
		while (seqIter.hasNext()) {
			RNASequence currentSequence = seqIter.next();
			String[] silvaTaxonomy = currentSequence.getTaxonomy().getSilvaTaxonomyPath().split(";");
			int result = findStringTax(taxSearch, silvaTaxonomy);
			if (result >= 0) {
				silvaCount++;
				result = -1;
			} else if (taxSearch.equalsIgnoreCase("Eukaryota")) {
				System.out.println(currentSequence.getTaxonomy().getSilvaTaxonomyPath());
			}
			String[] emblTaxonomy = currentSequence.getTaxonomy().getEMBLTaxonomyPath().split(";");
			result = findStringTax(taxSearch, emblTaxonomy);
			if (result >= 0) {
				emblCount++;
				result = -1;
			}
		}
		System.out.println("There are " + silvaCount + " occurrencies of " + taxSearch + " in the silva taxonomy.");
		System.out.println("There are " + emblCount + " occurrencies of " + taxSearch + " in the EMBL taxonomy.");
	}
	
	public static int findStringTax(String currentTax, String[] currentTaxonomy) {
		for (int i = currentTaxonomy.length-1; i >= 0; i--) {
			if (currentTaxonomy[i].equalsIgnoreCase(currentTax))
				return i;
		}
		return -1;
	}
	
	public void printStatisticsForAllRanks() {
		for (String currRank : rankFiles) {
			rankFile = currRank;
			currentRank = rankFile.substring(rankFile.lastIndexOf("_")+1);
			boolean allOk = readPhylaFiles(rankFile, false);
			if (allOk) {
				System.out.println("Ranks are read into the datastructure!");
				int begin = rankFile.lastIndexOf("_");
				String currentRank = rankFile.substring(begin);
				String phylaFile = phyloXMLFile + "_" + currentRank;
				System.out.println("Printing " + currentRank + " statistics...");
				printPhylaStatistics(phylaFile, false);
			}
			else
				System.out.println("The files were not read!");
		}
	}
	
	public void checkNonRedundant() {
		try {
			Scanner scanAmplicons = new Scanner(new FileInputStream(ampliconsFile));
			LinkedList<String> listAmplicons = new LinkedList<String>();
			while (scanAmplicons.hasNextLine()) {
				listAmplicons.add(scanAmplicons.nextLine());
			}
			Collections.sort(listAmplicons);
			ListIterator<RNASequence> nonRedIter = nonRedundantSequences.listIterator();
			LinkedList<String> nonRedAmplicons = new LinkedList<String>();
			while (nonRedIter.hasNext()) {
				RNASequence currentSequence = nonRedIter.next();
				nonRedAmplicons.add(currentSequence.getAccession().getID());
			}
			Collections.sort(nonRedAmplicons);
			
			PrintWriter amplicFile = new PrintWriter(ampliconsFile + "_originale");
			ListIterator<String> amplicIter = listAmplicons.listIterator();
			while (amplicIter.hasNext()) {
				amplicFile.println(amplicIter.next());
			}
			amplicFile.close();
			
			amplicFile = new PrintWriter(ampliconsFile + "_nonredun");
			ListIterator<String> nonRedAmplicIter = nonRedAmplicons.listIterator();
			while (nonRedAmplicIter.hasNext()) {
				amplicFile.println(nonRedAmplicIter.next());
			}
			amplicFile.close();
			scanAmplicons.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void readXMLFile(String xmlFile) {
		try {
			DocumentBuilderFactory factory = DocumentBuilderFactory.newInstance();
			factory.setIgnoringElementContentWhitespace(true);
			DocumentBuilder builder = factory.newDocumentBuilder();
			File file = new File(xmlFile);
			Document doc = builder.parse(file);
			Element root = doc.getDocumentElement();
			NodeList children = root.getChildNodes();
			for (int i = 0; i < children.getLength(); i++) {
				Node child = children.item(i);
				if (child instanceof Element) {
					Element childElement = (Element) child;
					NodeList grandChildren = childElement.getChildNodes();
					String sequenceData = null;
					String taxonomyIDString = null;
					String accessionIDString = null;
					String emblPath = null;
					String ncbiID = null;
					String silvaPath = null;
					String silvaName = childElement.getAttribute("name");
					for (int j = 0; j < grandChildren.getLength(); j++) {
						Node grandChild = grandChildren.item(j);
						if (grandChild instanceof Element) {
							Element grandChildElement = (Element) grandChild;
							if (grandChildElement.getTagName().equalsIgnoreCase("acc"))
								accessionIDString = grandChildElement.getTextContent();
							else if (grandChildElement.getTagName().equalsIgnoreCase("ALIGNMENT")) {
								Node grandGrandChild = grandChildElement.getFirstChild();
								if (grandGrandChild instanceof Element) {
									Element grandGrandChildElement = (Element) grandGrandChild;
									if (grandGrandChildElement.getTagName().equalsIgnoreCase("data"))
										sequenceData = grandGrandChildElement.getTextContent();
								}	
							}
							else if (grandChildElement.getTagName().equalsIgnoreCase("tax_embl"))
								emblPath = grandChildElement.getTextContent();
							else if (grandChildElement.getTagName().equalsIgnoreCase("tax_embl_name"))
								taxonomyIDString = grandChildElement.getTextContent();
							else if (grandChildElement.getTagName().equalsIgnoreCase("tax_xref_embl"))
								ncbiID = grandChildElement.getTextContent();
							else if (grandChildElement.getTagName().equalsIgnoreCase("tax_slv"))
								silvaPath = grandChildElement.getTextContent();
						}
					}
					RNASequence sequence = new RNASequence(sequenceData, AmbiguityRNACompoundSet.getRNACompoundSet());
					TaxonomyID taxonomyID = new TaxonomyID(taxonomyIDString, DataSource.ENA);
					taxonomyID.setEMBLTaxonomyPath(emblPath);
					taxonomyID.setSilvaTaxonomyPath(silvaPath);
					taxonomyID.setSilvaName(silvaName);
					AccessionID accessionID = new AccessionID(accessionIDString);
					accessionID.setNCBIId(ncbiID);
					sequence.setTaxonomy(taxonomyID);
					sequence.setAccession(accessionID);
					xmlSequences.add(sequence);
				}
			}
		} catch (Exception e) {
			e.printStackTrace();
		}
	}

	public int searchReversePrimerPosition(String rPrimer, String searchSequenceString, boolean mismatchAvoid) {
		RNASequence reversePrimerSequence = new RNASequence(rPrimer);
		String reverseIt = reversePrimerSequence.getInverse().getSequenceAsString();
		String stringPrimerSequence = new RNASequence(reverseIt).getComplement().getSequenceAsString();
		RNASequence primerSequence = new RNASequence(stringPrimerSequence);
		NucleotideCompound[] searchSequence = new NucleotideCompound[searchSequenceString.length()]; 
		(new ArrayListSequenceReader<NucleotideCompound>(searchSequenceString, AmbiguityRNACompoundSet.getRNACompoundSet())).getAsList().toArray(searchSequence);

		Iterator<NucleotideCompound> iterPrim = primerSequence.iterator();
		NucleotideCompound currentPrimerBase = null;
		NucleotideCompound currentSeqBase = null;
		for (int firstCounter = 0; firstCounter < searchSequence.length; firstCounter++) {
			currentSeqBase = searchSequence[firstCounter];
			currentPrimerBase = iterPrim.next();
			boolean localMismatchAvoid = mismatchAvoid;
			if (currentSeqBase.equalsIgnoreCase(currentPrimerBase) || !localMismatchAvoid) {
				int secondCounter = firstCounter;
				int thirdCounter = firstCounter;
				boolean found = false;
				while (secondCounter < searchSequence.length) {
					if (!(searchSequence[secondCounter].getBase().equals("-")) && !(searchSequence[secondCounter].getBase().equals(".")) ) {
						if (searchSequence[secondCounter].equalsIgnoreCase(currentPrimerBase) || !localMismatchAvoid) {
							if (iterPrim.hasNext()) {
								if (!searchSequence[secondCounter].equalsIgnoreCase(currentPrimerBase)) {
									localMismatchAvoid = true;
								}
								secondCounter++;
								currentPrimerBase = iterPrim.next();
								thirdCounter++;
							}
							else {
								secondCounter++;
								thirdCounter++;
								found = true;
							}
						}
						else {
							secondCounter = searchSequence.length;
						}
					}
					else {
						secondCounter++;
						thirdCounter++;
					}
					if (found)
						return thirdCounter;
				}
			}
			iterPrim = primerSequence.iterator();
		}
		return -1;
	}
	
	public int searchForwardPrimerPosition(String fPrimer, String searchSequenceString, boolean mismatchAvoid) {
		RNASequence primerSequence = new RNASequence(fPrimer);
		NucleotideCompound[] searchSequence = new NucleotideCompound[searchSequenceString.length()]; 
		(new ArrayListSequenceReader<NucleotideCompound>(searchSequenceString, AmbiguityRNACompoundSet.getRNACompoundSet())).getAsList().toArray(searchSequence);		
		Iterator<NucleotideCompound> iterPrim = primerSequence.iterator();
		NucleotideCompound currentPrimerBase = null;
		NucleotideCompound currentSeqBase = null;
		for (int firstCounter = 0; firstCounter < searchSequence.length; firstCounter++) {
			currentSeqBase = searchSequence[firstCounter];
			currentPrimerBase = iterPrim.next();
			boolean localMismatchAvoid = mismatchAvoid;
			if (currentSeqBase.equalsIgnoreCase(currentPrimerBase) || !localMismatchAvoid) {
				int secondCounter = firstCounter;
				while (secondCounter < searchSequence.length) {
					if (!(searchSequence[secondCounter].getBase().equals("-")) && !(searchSequence[secondCounter].getBase().equals(".")) ) {
						if (searchSequence[secondCounter].equalsIgnoreCase(currentPrimerBase) || !localMismatchAvoid) {
							if (!searchSequence[secondCounter].equalsIgnoreCase(currentPrimerBase)) {
								localMismatchAvoid = true;
							}
							secondCounter++;
							currentPrimerBase = iterPrim.next();
						}
						else {
							secondCounter = searchSequence.length;
						}
					}
					else {
						secondCounter++;
					}
					if (!iterPrim.hasNext())
						return firstCounter;
				}
			}
			iterPrim = primerSequence.iterator();
		}
		return -1;
	}
	
	public int generateAmplicons(boolean mismatchAvoid) {
		ListIterator<RNASequence> sequencesIterator = xmlSequences.listIterator();
		int numberAmplicons = 0;
		while (sequencesIterator.hasNext()) {
			RNASequence currentSequence = sequencesIterator.next();
			String currentSequenceString = currentSequence.getSequenceAsString();
			int indexFor = searchForwardPrimerPosition(forwardPrimer, currentSequenceString, mismatchAvoid);
			int indexRev = 0;
			if (indexFor > 0)
				indexRev = searchReversePrimerPosition(reversePrimer, currentSequenceString, mismatchAvoid);
			if (indexRev > 0 && indexRev > indexFor) {
				numberAmplicons++;
				System.out.println("Current number of amplicons: " + numberAmplicons);
				String amplicon = currentSequenceString.substring(indexFor, indexRev);
				String removedGaps = amplicon.replaceAll("-", "");
				
				RNASequence ampliconRNA = new RNASequence(removedGaps, AmbiguityRNACompoundSet.getRNACompoundSet());
				ampliconRNA.setTaxonomy(currentSequence.getTaxonomy());
				ampliconRNA.setAccession(currentSequence.getAccession());
				ampliconSequencesRNA.add(ampliconRNA);
				
				String ampliconDNAString = removedGaps.replaceAll("U", "T");
				
				DNASequence ampliconDNA1 = new DNASequence(ampliconDNAString, AmbiguityDNACompoundSet.getDNACompoundSet());
				ampliconDNA1.setTaxonomy(currentSequence.getTaxonomy());
				String accessionIDString = currentSequence.getAccession().getID();
				accessionIDString += "+";
				AccessionID accessionID1 = new AccessionID(accessionIDString);
				accessionID1.setNCBIId(currentSequence.getAccession().getNCBIId());
				ampliconDNA1.setAccession(accessionID1);
				ampliconSequences.add(ampliconDNA1);
				
				Iterator<NucleotideCompound> seqIter = ampliconDNA1.iterator();
				String complSequenceString = "";
				while (seqIter.hasNext()) {
					NucleotideCompound currentCompound = seqIter.next();
					complSequenceString += currentCompound.getComplement().getShortName();
				}
				DNASequence ampliconDNA2 = new DNASequence(complSequenceString, AmbiguityDNACompoundSet.getDNACompoundSet());
				ampliconDNA2.setTaxonomy(ampliconDNA1.getTaxonomy());
				accessionIDString = currentSequence.getAccession().getID();
				accessionIDString += "-";
				AccessionID accessionID2 = new AccessionID(accessionIDString);
				accessionID2.setNCBIId(accessionID1.getNCBIId());
				ampliconDNA2.setAccession(accessionID2);
				ampliconSequences.add(ampliconDNA2);
			}
		}
		return numberAmplicons;
	}
	
	public void generateAllSequencesInSpecificRegion(String markString, int start, int stop) {
		PrintWriter fullLengthSeqWriter = null;
		PrintWriter lengthWriter = null;
		try {
			fullLengthSeqWriter = new PrintWriter(xmlFile.substring(0,xmlFile.lastIndexOf(".")) + markString + ".fasta");
			lengthWriter = new PrintWriter(xmlFile.substring(0, xmlFile.lastIndexOf(".")) + markString + "_length");
			ListIterator<RNASequence> seqIter = xmlSequences.listIterator();
			while (seqIter.hasNext()) {
				RNASequence currentSequence =seqIter.next();
				String currentWholeSeqString = currentSequence.getSequenceAsString();
				if (currentWholeSeqString.length() == 24782) {
					String current1800String = currentWholeSeqString.substring(start, stop);
					String removedEndGaps = current1800String.replaceAll("\\.", "");
					String removedGaps = removedEndGaps.replaceAll("-", "");
					String printMe = ">";
					printMe += currentSequence.getAccession().getID();
					fullLengthSeqWriter.println(printMe);
					fullLengthSeqWriter.println(removedGaps);
					lengthWriter.println(removedGaps.length());
				} else {
					System.out.println("The sequences are not 24782 nucleotides long but " + currentWholeSeqString.length());
					System.exit(0);
				}
			}
			fullLengthSeqWriter.close();
			lengthWriter.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void printAllToFasta() {
		PrintWriter allFasta = null;
		try {
			allFasta = new PrintWriter(xmlFile.substring(0, xmlFile.lastIndexOf(".")) + ".fasta");
			ListIterator<RNASequence> seqIter = xmlSequences.listIterator();
			while (seqIter.hasNext()) {
				RNASequence currentSequence = seqIter.next();
				String currentSequenceString = currentSequence.getSequenceAsString();
				String removedEndGaps = currentSequenceString.replaceAll("\\.", "");
				String removedGaps = removedEndGaps.replaceAll("-", "");
				String printMe = ">";
				printMe += currentSequence.getAccession().getID();
				allFasta.println(printMe);
				allFasta.println(removedGaps);
			}
			allFasta.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void printAmplicons(String amplicFile, SequenceType seqType, boolean printSense) {
		PrintWriter amplic = null;
		PrintWriter amplicLength = null;
		PrintWriter amplicNames = null;
		try {
			amplic = new PrintWriter(amplicFile);
			if (seqType.equals(SequenceType.RNA)) {
				amplicLength = new PrintWriter(fileAmpliconsLength);
				amplicNames = new PrintWriter(fileAmpliconsNames);
				ListIterator<RNASequence> sequencesIterator = ampliconSequencesRNA.listIterator();
				while (sequencesIterator.hasNext()) {
					RNASequence currentSequence = sequencesIterator.next();
					String currentSequenceString = currentSequence.getSequenceAsString();
					String printMe = ">";
					printMe += currentSequence.getAccession().getID();
					amplic.println(printMe);
					amplic.println(currentSequenceString);
					amplicLength.println(currentSequenceString.length());
					amplicNames.println(currentSequence.getAccession().getID());
				}
				amplicLength.close();
				amplicNames.close();
			}
			else if (seqType.equals(SequenceType.DNA)){
				ListIterator<DNASequence> sequencesIterator = ampliconSequences.listIterator();
				while (sequencesIterator.hasNext()) {
					DNASequence currentSequence = sequencesIterator.next();
					String currentSequenceString = currentSequence.getSequenceAsString();
					String printMe = ">";
					String name = null;
					if (!printSense) 
						name = currentSequence.getAccession().getID().substring(0, currentSequence.getAccession().getID().length()-1);
					else
						name = currentSequence.getAccession().getID();
					printMe += name;
					amplic.println(printMe);
					amplic.println(currentSequenceString);
				}
			}
			else {
				System.out.println("There is no such Sequence Type!");
				System.exit(1);
			}
			amplic.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void printSilva2NCBIMap(String printFile) {
		PrintWriter silvaNCBI = null;
		try {
			silvaNCBI = new PrintWriter(printFile);
			ListIterator<RNASequence> sequencesIterator = xmlSequences.listIterator();
			while (sequencesIterator.hasNext()) {
				RNASequence currentSequence = sequencesIterator.next();
				silvaNCBI.println(currentSequence.getAccession().getID() + "\t" + currentSequence.getAccession().getNCBIId());
			}
			silvaNCBI.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void constructSequencesSites() {
		ListIterator<RNASequence> sequencesIterator = xmlSequences.listIterator();
		while (sequencesIterator.hasNext()) {
			RNASequence currentSequence = sequencesIterator.next();
			if (currentSequence != null) {
				Iterator<NucleotideCompound> currSeqIter =currentSequence.iterator();
				if (allSequenceSites.isEmpty()) {
					while (currSeqIter.hasNext()) {
						currSeqIter.next();
						allSequenceSites.add(new Site());
					}
					currSeqIter =currentSequence.iterator();
				}
				Iterator<Site> siteIter = allSequenceSites.iterator();
				while (currSeqIter.hasNext()) {
					NucleotideCompound currentCompound = currSeqIter.next();
					Site currentSite = siteIter.next();
					if (currentSite.getLastOne() != null) {
						if (!currentCompound.equalsIgnoreCase(currentSite.getLastOne())) {
							currentSite.increaseChange();
						}
					}
					currentSite.setlastOne(currentCompound);
					Integer currentInteger = currentSite.getIntegerValue(currentCompound);
					if (currentInteger == null)
						currentInteger = 0;
					currentInteger++;
					currentSite.addCompound(currentCompound, currentInteger);
				}
			}
		}
	}
	
	public RNASequence searchSequenceBySilvaName() {
		ListIterator<RNASequence> sequencesIterator = xmlSequences.listIterator();
		while (sequencesIterator.hasNext()) {
			RNASequence currentSequence = sequencesIterator.next();
			if (currentSequence.getTaxonomy().getSilvaName().equalsIgnoreCase(referenceSilvaName))
				return currentSequence;
		}
		return null;
	}
	
	public void buildNonRedundantSequences(String ampliconsFile) {
		try {
			nonRedundantSequences.clear();
			Scanner scan = new Scanner(new FileInputStream(ampliconsFile));
			LinkedList<String> listAmplicons = new LinkedList<String>();
			while (scan.hasNextLine()) {
				listAmplicons.add(scan.nextLine());
			}
			ListIterator<RNASequence> sequencesIterator = xmlSequences.listIterator();
			while (sequencesIterator.hasNext()) {
				RNASequence currentSequence = sequencesIterator.next();
				ListIterator<String> ampliconsIterator = listAmplicons.listIterator();
				boolean found = false;
				while (ampliconsIterator.hasNext() && !found) {
					if (currentSequence.getAccession().getID().equalsIgnoreCase(ampliconsIterator.next())) {
						nonRedundantSequences.add(currentSequence);
						found = true;
					}
				}
			}
			scan.close();
			System.out.println("There are " + listAmplicons.size() + " amplicons in the list.");
			System.out.println(nonRedundantSequences.size() + " of them are found in the database.");
			System.out.println("There are " + xmlSequences.size() + " sequences in the database.");
		} catch (FileNotFoundException fException) {
			fException.printStackTrace();
		}
	}

	public void readPhyloXMLFile(String phyloXMLFile) {
		try {
			PhyloXmlParser treeParser = new PhyloXmlParser();
			treeParser.setSource(new File(phyloXMLFile));
			phylogenies = treeParser.parse();
		} catch (IOException e) {
			System.out.println("File " + phyloXMLFile + " not found!");
			e.printStackTrace();
		}
	}
	
	public boolean readPhylaFiles(String rankFile, boolean sub) {
		boolean filesRead = false;
		try {
			rankList.clear();
			classifiedRanks.clear();
			subphylaList.clear();
			silvaRankList.clear();
			chosenRankList.clear();
			lowerRankList.clear();
			Scanner scan = new Scanner(new FileInputStream(rankFile));
			while (scan.hasNextLine()) {
				String phylaString = scan.nextLine();
				rankList.add(phylaString);
			}
			scan.close();
			if (sub) {
				scan = new Scanner(new FileInputStream("subphyla"));
				while (scan.hasNextLine()) {
					String line = scan.nextLine();
					String[] pairs = line.split(" ");
					if (pairs.length < 2)
						System.exit(0);
					subphylaList.put(pairs[0], pairs[1]);
				}
				scan.close();
			}
			if (!rankList.isEmpty() || !subphylaList.isEmpty())
				filesRead = true;
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
		return filesRead;
	}
	
	public void printPhylaStatistics(String phylaFile, boolean sub) {
		PhylogenyNodeIterator phylogenyIter = phylogenies[0].iteratorPostorder();
		chosenRankList.clear();
		if (!chosenRankList.isEmpty())
			System.exit(0);
		LinkedList<Phylum> chosenSubPhylaList = new LinkedList<Phylum>();
		try {
			FileWriter outFile = new FileWriter(phylaFile);
			PrintWriter out = new PrintWriter(outFile);
			while (phylogenyIter.hasNext()) {
				PhylogenyNode currentNode = phylogenyIter.next();
				ListIterator<String> phylaIter = rankList.listIterator();
				boolean found = false;
				while (phylaIter.hasNext() && !found) {
					String currentPhylum = phylaIter.next();
					if (!currentNode.isExternal() && currentNode.getNodeData().getTaxonomy().getScientificName().equalsIgnoreCase(currentPhylum)) {
						ListIterator<Phylum> chosenPhylaIter = chosenRankList.listIterator();
						boolean foundInChosen = false;
						while (chosenPhylaIter.hasNext() && !foundInChosen) {
							int currentIndex = chosenPhylaIter.nextIndex();
							Phylum currentChosenPhylum = chosenPhylaIter.next();
							if (currentChosenPhylum.getName().equalsIgnoreCase(currentPhylum)) {
								foundInChosen = true;
								if (currentChosenPhylum.getDepth() > PhylogenyMethods.calculateDepth(currentNode)) {
									chosenRankList.set(currentIndex, new Phylum(currentPhylum, PhylogenyMethods.calculateDepth(currentNode), 
											currentNode.getNumberOfExternalNodes(), PhylogenyMethods.calculateSumOfDistinctTaxonomies(currentNode),
											currentNode));
								}
							}
						}
						if (!foundInChosen) {
							chosenRankList.add(new Phylum(currentPhylum, PhylogenyMethods.calculateDepth(currentNode), currentNode.getNumberOfExternalNodes(), 
								PhylogenyMethods.calculateSumOfDistinctTaxonomies(currentNode), currentNode));
						}
						found = true;
					}
				}
			}
			if (sub) {
				phylogenyIter = phylogenies[0].iteratorPostorder();
				while (phylogenyIter.hasNext()) {
					PhylogenyNode currentNode = phylogenyIter.next();
					if (subphylaList.containsKey(currentNode.getNodeData().getTaxonomy().getScientificName())) {
						ListIterator<Phylum> chosenSubPhylaIter = chosenSubPhylaList.listIterator();
						boolean foundInChosen = false;
						while (chosenSubPhylaIter.hasNext() && !foundInChosen) {
							int currentIndex = chosenSubPhylaIter.nextIndex();
							Phylum currentChosenSubPhylum = chosenSubPhylaIter.next();
							if (currentChosenSubPhylum.getName().equalsIgnoreCase(currentNode.getNodeData().getTaxonomy().getScientificName())) {
								foundInChosen = true;
								if (currentChosenSubPhylum.getDepth() > PhylogenyMethods.calculateDepth(currentNode)) {
									chosenSubPhylaList.set(currentIndex, new Phylum(currentNode.getNodeData().getTaxonomy().getScientificName(),
											PhylogenyMethods.calculateDepth(currentNode), currentNode.getNumberOfExternalNodes(),
											PhylogenyMethods.calculateSumOfDistinctTaxonomies(currentNode), currentNode));
								}
							}
						}
						if (!foundInChosen) {
							chosenSubPhylaList.add(new Phylum(currentNode.getNodeData().getTaxonomy().getScientificName(),
									PhylogenyMethods.calculateDepth(currentNode), currentNode.getNumberOfExternalNodes(),
									PhylogenyMethods.calculateSumOfDistinctTaxonomies(currentNode), currentNode));
						}
					}
				}
				
				ListIterator<Phylum> chosenSubPhylaIter = chosenSubPhylaList.listIterator();
				while (chosenSubPhylaIter.hasNext()) {
					Phylum currentSubPhylum = chosenSubPhylaIter.next();
					if (subphylaList.containsKey(currentSubPhylum.getName())) {
						ListIterator<Phylum> chosenPhylaIter = chosenRankList.listIterator();
						boolean exists = false;
						while (chosenPhylaIter.hasNext() && !exists) {
							Phylum currentChosenPhylum = chosenPhylaIter.next();
							if (currentChosenPhylum.getName().equalsIgnoreCase(subphylaList.get(currentSubPhylum.getName()))) {
								exists = true;
								int currentSumOfTax = currentChosenPhylum.getSumOfDistinctTaxonomies();
								currentChosenPhylum.setSumOfDistinctTaxonomies(currentSumOfTax + currentSubPhylum.getSumOfDistinctTaxonomies());
								currentChosenPhylum.addSubCladeSource(currentSubPhylum);
							}
						}
						if (!exists) {
							chosenRankList.add(new Phylum(subphylaList.get(currentSubPhylum.getName()), currentSubPhylum.getDepth(), 
									currentSubPhylum.getNumberOfExternalNodes(), currentSubPhylum.getSumOfDistinctTaxonomies(), currentSubPhylum,
									currentSubPhylum.getPhylogenyNode()));
						}
					}
				}
			}
			
			Collections.sort(chosenRankList);
			ListIterator<Phylum> chosenPhylaIter = chosenRankList.listIterator();
			while(chosenPhylaIter.hasNext()) {
				Phylum chosenPhylum = chosenPhylaIter.next();
				out.println(chosenPhylum.getName() + "\t\t\t" + chosenPhylum.getDepth() + "\t" + chosenPhylum.getNumberOfExternalNodes() + "\t" 
						+ chosenPhylum.getSumOfDistinctTaxonomies());
				if (sub && !chosenPhylum.getSubClades().isEmpty()) {
					ListIterator<Phylum> subCladesIter = chosenPhylum.getSubClades().listIterator();
					while (subCladesIter.hasNext()) {
						out.println("\t" + subCladesIter.next().getName());
					}
				}
			}
			out.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}
	
	public void createAndPrintAllPhylumAmplicons() {
		ListIterator<String> allPhylaIter = rankList.listIterator();
		while (allPhylaIter.hasNext()) {
			long startPhylum = System.currentTimeMillis();
			searchPhylum = allPhylaIter.next();
			Phylogeny[] subtrees = null;
			phylumAmplicons.clear();
			ListIterator<Phylum> chosenPhylaIter = chosenRankList.listIterator();
			Phylum currentPhylum = null;
			boolean found = false;
			while (chosenPhylaIter.hasNext() && !found) {
				currentPhylum = chosenPhylaIter.next();
				if (currentPhylum.getName().equalsIgnoreCase(searchPhylum)) {
					found = true;
				}
			}
			if (found) {
				subtrees = new Phylogeny[currentPhylum.getSubClades().size()+1];
				subtrees[0] = phylogenies[0].subTree(currentPhylum.getPhylogenyNode());
				ListIterator<Phylum> subCladesIter = currentPhylum.getSubClades().listIterator();
				int subtreeCounter = 1;
				while (subCladesIter.hasNext()) {
					Phylum currSubPhylum = subCladesIter.next();
					subtrees[subtreeCounter] = phylogenies[0].subTree(currSubPhylum.getPhylogenyNode());
					subtreeCounter++;
				}
								
				assert(subtrees.length == (currentPhylum.getSubClades().size() + 1));
				ListIterator<RNASequence> amplicIter = ampliconSequencesRNA.listIterator();
				while (amplicIter.hasNext()) {
					RNASequence currentAmplicon = amplicIter.next();
					PhylogenyNodeIterator phylogenyIter = subtrees[0].iteratorExternalForward();
					boolean foundAmplicon = false;
					while (phylogenyIter.hasNext() && !foundAmplicon) {
						PhylogenyNode currentNode = phylogenyIter.next();
						if (currentNode.getNodeData().getTaxonomy().getScientificName().equalsIgnoreCase(currentAmplicon.getTaxonomy().getID()) &&
								currentNode.getNodeData().getSequence().getAccession().getValue().equalsIgnoreCase(currentAmplicon.getAccession().getID())) {
							phylumAmplicons.add(currentAmplicon);
							foundAmplicon = true;
						}
					}
				}
				PrintWriter amplic = null;
				PrintWriter amplicLength = null;
				String amplicFile = ampliconsFile + "_" + searchPhylum;
				try {
					amplic = new PrintWriter(amplicFile);
					amplicLength = new PrintWriter(amplicFile + "_length");
					ListIterator<RNASequence> sequencesIterator = phylumAmplicons.listIterator();
					@SuppressWarnings("unused")
					int sequencesCounter = 0;
					while (sequencesIterator.hasNext()) {
						RNASequence currentSequence = sequencesIterator.next();
						sequencesCounter++;
						String currentSequenceString = currentSequence.getSequenceAsString();
						String printMe = ">";
						printMe += currentSequence.getAccession().getID();
						amplic.println(printMe);
						amplic.println(currentSequenceString);
						amplicLength.println(currentSequenceString.length());
					}
					amplic.close();
					amplicLength.close();
				} catch (FileNotFoundException e) {
					e.printStackTrace();
				}
			}
			else {
				PrintWriter amplic = null;
				String amplicFile = ampliconsFile + "_" + searchPhylum + "_absent";
				try {
					amplic = new PrintWriter(amplicFile);
					amplic.println("This phylum is not present!");
					amplic.close();
				} catch (FileNotFoundException e) {
					e.printStackTrace();
				}
			}
			long endPhylum = System.currentTimeMillis();
			System.out.println("Phylum " + searchPhylum + " printed to file! \t" + (double)(endPhylum-startPhylum)/1000 + " s");
		}
		phylumAmplicons.clear();
	}
	
//	public void createSubtreeAndPrintPhylumStatistics(String searchPhylum) {
		
//	}
	
	public void findTaxonomy() {
		ListIterator<RNASequence> seqIterator = nonRedundantSequences.listIterator();
		try {
			PrintWriter taxWriter = new PrintWriter(ampliconsFile + "_taxonomies");
			while (seqIterator.hasNext()) {
				RNASequence currentSequence = seqIterator.next();
				taxWriter.println(">" + currentSequence.getAccession().getID());
				taxWriter.println(currentSequence.getTaxonomy().getSilvaTaxonomyPath());
			}
			taxWriter.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
/*	public String searchTaxonomy(String searchQuery) {
		ListIterator<RNASequence> seqIter = xmlSequences.listIterator();
		String taxonomy = "";
		//String ncbiAccession = "";
		boolean found = false;
		while (seqIter.hasNext() && !found) {
			RNASequence currentSequence = seqIter.next();
			if (currentSequence.getAccession().getID().equalsIgnoreCase(searchQuery)) {
				taxonomy = currentSequence.getTaxonomy().getEMBLTaxonomyPath();
				//ncbiAccession = currentSequence.getAccession().getNCBIId();
				found = true;
			}
		}
		if (!found) {
			//System.out.println(searchQuery + " not found!");
			return null;
		}
		else {
			//System.out.println("The taxonomy of " + searchQuery + " is: " + taxonomy);
			//System.out.println("NCBI accession number is: " + ncbiAccession);
			return taxonomy;
		}
	}
*/	
	public void readNCBITaxonomy() {
		ncbiTaxonomies = new NCBITaxonomy[1092595];
		try {
			Scanner scanNames = new Scanner(new FileInputStream("names.dmp"));
			Scanner scanNodes = new Scanner(new FileInputStream("nodes.dmp"));
			while (scanNames.hasNextLine()) {
				String currentLine = scanNames.nextLine();
				String[] lineFields = currentLine.split("\\|");
				Integer tax_id = new Integer(lineFields[0].replaceAll("\t", ""));
				String name_txt = lineFields[1].replaceAll("\t", "");
				//String unique_name = lineFields[2].replaceAll("\t", "");
				String name_class = lineFields[3].replaceAll("\t", "");
				if (name_class.equalsIgnoreCase("scientific name")) {
					ncbiTaxonomies[tax_id] = new NCBITaxonomy(tax_id, name_txt);
					String currentNodeLine = scanNodes.nextLine();
					String[] lineFieldsNode = currentNodeLine.split("\\|");
					Integer tax_id_node = new Integer(lineFieldsNode[0].replaceAll("\t", ""));
					if (tax_id_node.equals(tax_id)) {
						Integer parent_tax_id_node = new Integer(lineFieldsNode[1].replaceAll("\t", ""));
						String rank = lineFieldsNode[2].replaceAll("\t", "");
						ncbiTaxonomies[tax_id].setParentUIDAndRank(parent_tax_id_node, rank);
					} else {
						System.out.println("Tax ids differ! " + tax_id + "!=" + tax_id_node);
						System.exit(0);
					}
				}
			}
			scanNames.close();
			scanNodes.close();
			Scanner scanMerged = new Scanner(new FileInputStream("merged.dmp"));
			while (scanMerged.hasNextLine()) {
				String currentLine = scanMerged.nextLine();
				String[] lineFields = currentLine.split("\\|");
				Integer tax_id_old = new Integer(lineFields[0].replaceAll("\t", ""));
				Integer tax_id_new = new Integer(lineFields[1].replaceAll("\t", ""));
				if (ncbiTaxonomies[tax_id_old] == null)
					ncbiTaxonomies[tax_id_old] = new NCBITaxonomy(tax_id_old);
				ncbiTaxonomies[tax_id_old].setInvalid(tax_id_new);
			}
			scanMerged.close();
			Scanner scanDeleted = new Scanner(new FileInputStream("delnodes.dmp"));
			while (scanDeleted.hasNextLine()) {
				String currentLine = scanDeleted.nextLine();
				String[] lineFields = currentLine.split("\\|");
				Integer tax_id_del = new Integer(lineFields[0].replaceAll("\t", ""));
				if (ncbiTaxonomies[tax_id_del] == null)
					ncbiTaxonomies[tax_id_del] = new NCBITaxonomy(tax_id_del);
				ncbiTaxonomies[tax_id_del].setDeleted();
			}
			scanDeleted.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void findRanksInNCBI(String rank, boolean all, int retainUID) {
		try {
			PrintWriter ncbiTaxWriter;
			LinkedList<NCBITaxonomy> chosenTax = new LinkedList<NCBITaxonomy>();
			LinkedList<NCBITaxonomy> UIDsToBeRemoved = new LinkedList<NCBITaxonomy>();
			for (NCBITaxonomy ncbiTax : ncbiTaxonomies) {
				if (ncbiTax != null && !ncbiTax.isInvalid() && !ncbiTax.isDeleted() && ncbiTax.getRank().equalsIgnoreCase(rank)) {
					chosenTax.add(ncbiTax);
				}
			}
			if (!all) {
				ListIterator<NCBITaxonomy> chosenTaxIter = chosenTax.listIterator();
				while (chosenTaxIter.hasNext()) {
					NCBITaxonomy currentTax = chosenTaxIter.next();
					NCBITaxonomy newCurrentTax = currentTax;
					boolean found = false;
					while (!found && newCurrentTax.getParentUID() != 1) {
						int parentUID = newCurrentTax.getParentUID();
						newCurrentTax = ncbiTaxonomies[parentUID];
						if (newCurrentTax.getUID() == retainUID) {
							found = true;
						}
					}
					if (!found) {
						UIDsToBeRemoved.add(currentTax);
					}
				}
			}
			if (all) {
				ncbiTaxWriter = new PrintWriter("all_ncbi_" + rank);
			}
			else {
				ListIterator<NCBITaxonomy> removeIter = UIDsToBeRemoved.listIterator();
				while (removeIter.hasNext()) {
					chosenTax.remove(removeIter.next());
				}
				ncbiTaxWriter = new PrintWriter(retainUID + "_ncbi_" + rank);
			}
			Collections.sort(chosenTax);
			ListIterator<NCBITaxonomy> chosenTaxIter = chosenTax.listIterator();
			while (chosenTaxIter.hasNext()) {
				NCBITaxonomy currentTax = chosenTaxIter.next();
				//ncbiTaxWriter.println(currentTax.getSciName() + "\t" + currentTax.getRank());
				ncbiTaxWriter.println(currentTax.getSciName());
			}
			ncbiTaxWriter.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void findRanksInSilva(String rank) {
		ListIterator<RNASequence> seqIter = xmlSequences.listIterator();
		while (seqIter.hasNext()) {
			RNASequence currentSequence = seqIter.next();
			LinkedList<String> validRanksOfSeq = new LinkedList<String>();
			Integer currentUID = new Integer(currentSequence.getAccession().getNCBIId());
			Integer seqID = currentUID;
			if (ncbiTaxonomies[currentUID] != null && ncbiTaxonomies[currentUID].isInvalid()) {
				currentUID = ncbiTaxonomies[currentUID].getMergedUID();
			}
			if (ncbiTaxonomies[currentUID] != null && !ncbiTaxonomies[currentUID].isInvalid() && !ncbiTaxonomies[currentUID].isDeleted()) {				
				Integer currentParentUID = ncbiTaxonomies[currentUID].getParentUID();
				String currentRank = ncbiTaxonomies[currentUID].getRank();
				if (!currentRank.equalsIgnoreCase("no rank"))
					validRanksOfSeq.add(currentRank);
				Integer currentHighestValidRankID = null;
			//	String currentSciName = ncbiTaxonomies[currentUID].getSciName();
			//	System.out.println(currentUID + "\t" + currentRank + "\t" + rank);
				if (!currentRank.equals(rank)) {
					boolean found = false;
					while (!found && currentParentUID != 1) {
						if (ncbiTaxonomies[currentParentUID] != null && !ncbiTaxonomies[currentParentUID].isInvalid() && !ncbiTaxonomies[currentUID].isDeleted()) {
							currentUID = ncbiTaxonomies[currentParentUID].getUID();
							currentRank = ncbiTaxonomies[currentUID].getRank();
							if (!currentRank.equalsIgnoreCase("no rank"))
								validRanksOfSeq.add(currentRank);
							if (currentRank.equalsIgnoreCase(rank)) {
								found = true;
							} else {
								if (!currentRank.equalsIgnoreCase("no rank") && currentHighestValidRankID == null)
									currentHighestValidRankID = currentUID;
								currentParentUID = ncbiTaxonomies[currentUID].getParentUID();
							}
						} else if (ncbiTaxonomies[currentParentUID].isInvalid()) {
							currentUID = ncbiTaxonomies[currentParentUID].getMergedUID();
							currentRank = ncbiTaxonomies[currentUID].getRank();
							if (!currentRank.equalsIgnoreCase("no rank"))
								validRanksOfSeq.add(currentRank);
							if (currentRank.equalsIgnoreCase(rank)) {
								found = true;
							} else {
								if (!currentRank.equalsIgnoreCase("no rank") && currentHighestValidRankID == null)
									currentHighestValidRankID = currentUID;
								currentParentUID = ncbiTaxonomies[currentUID].getParentUID();
							}
						}
					}
					if (found && !silvaRankList.contains(ncbiTaxonomies[currentUID].getSciName())) {
						silvaRankList.add(ncbiTaxonomies[currentUID].getSciName());
					} else if (!found && !lowerRankList.containsKey(ncbiTaxonomies[currentHighestValidRankID].getSciName())){
						lowerRankList.put(seqID, validRanksOfSeq);
					}
				} else if (!silvaRankList.contains(ncbiTaxonomies[currentUID].getSciName())) {
					silvaRankList.add(ncbiTaxonomies[currentUID].getSciName());
				}
			}
		}
		Collections.sort(silvaRankList);
		try {
			PrintWriter rankWrite = new PrintWriter("all_silva_" + rank);
			ListIterator<String> rankIter = silvaRankList.listIterator();
			while (rankIter.hasNext()) {
				rankWrite.println(rankIter.next());
			}
			rankWrite.close();
			rankWrite = new PrintWriter("all_sub_silva_" + rank);
			for (Map.Entry<Integer, LinkedList<String>> entry : lowerRankList.entrySet()) {
				rankWrite.println(entry.getKey() + "\t" + entry.getValue().toString());
			}
			rankWrite.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void readClassifiedFasta(String inputFasta, String currRank) {
		try {
			classifiedRanks.clear();
			PrintWriter notfoundWriter = new PrintWriter(inputFasta.substring(0, inputFasta.indexOf("_")) + "_" + currRank + "_notfound");
			Scanner scanClassResults = new Scanner(new FileInputStream(inputFasta));
			while (scanClassResults.hasNextLine()) {
				String currentLine = scanClassResults.nextLine();
				if (currentLine.startsWith(">")) {
					String[] header = currentLine.split("\t");
					String[] taxonomy = header[1].split(";");
					ListIterator<String> rankIter = rankList.listIterator();
					boolean found = false;
					String currentRank = null;
					while (!found && rankIter.hasNext()) {
						currentRank = rankIter.next();
						for (String classTax : taxonomy) {
							if (classTax.equals(currentRank))
								found = true;
						}
					}
					if (found) {
						if (classifiedRanks.containsKey(currentRank)) {
							int currentValue = classifiedRanks.get(currentRank);
							classifiedRanks.put(currentRank, currentValue+1);
						}
						else {
							classifiedRanks.put(currentRank, 1);
						}
					}
					else {
						for (int taxCounter = 2; taxCounter < taxonomy.length; taxCounter++) {
							notfoundWriter.print(taxonomy[taxCounter] + ";");
						}
						notfoundWriter.println();
						if (classifiedRanks.containsKey("not found")) {
							int currentValue = classifiedRanks.get("not found");
							classifiedRanks.put("not found", currentValue+1);
						}
						else {
							classifiedRanks.put("not found", 1);
						}
					}
				}
			}
			scanClassResults.close();
			notfoundWriter.close();
			TreeMap<String, Integer> sortedRanks = new TreeMap<String, Integer>(classifiedRanks);
			PrintWriter sortedRanksWriter = new PrintWriter(inputFasta.substring(0, inputFasta.indexOf("_")) + "_" + currRank);
			String firstRank = sortedRanks.firstKey();
			if (!firstRank.equals("not found"))
				sortedRanksWriter.println(firstRank + "\t" + sortedRanks.get(firstRank));
			while (!firstRank.equals(sortedRanks.lastKey())) {
				firstRank = sortedRanks.higherKey(firstRank);
				if (!firstRank.equals("not found"))
					sortedRanksWriter.println(firstRank + "\t" + sortedRanks.get(firstRank));
			}
			if (sortedRanks.containsKey("not found"))
				sortedRanksWriter.println("not found" + "\t" + sortedRanks.get("not found"));
			sortedRanksWriter.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void compareTwoSilvaDatabases(String originale, String reformatio) {
		try {
			Scanner scanTodo = new Scanner(new FileInputStream(originale));
			String accessionID = null;
			while (scanTodo.hasNextLine()) {
				String currentLine = scanTodo.nextLine();
				if (currentLine.startsWith(">")) {
					accessionID = currentLine.substring(1);
					databaseOriginale.add(accessionID);
				}
			}
			Collections.sort(databaseOriginale);
			scanTodo.close();
			scanTodo = new Scanner(new FileInputStream(reformatio));
			accessionID =null;
			while (scanTodo.hasNextLine()) {
				String currentLine = scanTodo.nextLine();
				if (currentLine.startsWith(">")) {
					accessionID = currentLine.substring(1);
					databaseReformatio.add(accessionID);
				}
			}
			Collections.sort(databaseReformatio);
			scanTodo.close();
			ListIterator<String> reformIter = databaseReformatio.listIterator();
			while (reformIter.hasNext()) {
				String currentReformString = reformIter.next();
				boolean found = false;
				ListIterator<String> origIter = databaseOriginale.listIterator();
				while (!found && origIter.hasNext()) {
					int currentOrigIndex = origIter.nextIndex();
					String currentOrigString = origIter.next();
					if (currentReformString.equalsIgnoreCase(currentOrigString)) {
						databaseOriginale.remove(currentOrigIndex);
						found = true;
					}
					//if (currentReformString.compareToIgnoreCase(currentOrigString) > 0)
					//	found = true;
				}
			}
			PrintWriter printRestOriginale = new PrintWriter("Database_originale");
			ListIterator<String> origRestIter = databaseOriginale.listIterator();
			while (origRestIter.hasNext()) {
				printRestOriginale.println(origRestIter.next());
			}
			printRestOriginale.close();
			PrintWriter printReformatio = new PrintWriter("Database_reformatio");
			ListIterator<String> reformatioIter = databaseReformatio.listIterator();
			while (reformatioIter.hasNext()) {
				printReformatio.println(reformatioIter.next());
			}
			printReformatio.close();
			System.out.println(databaseOriginale.size() + " sequences have been removed from the original database");
			
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
	public void printNamesAndTaxa() {
		PrintWriter nameRedund = null;
		PrintWriter taxaRedund = null;
		try {
			LinkedList<String> redundNames = new LinkedList<String>();
			LinkedList<String> redundTaxa = new LinkedList<String>();
			ListIterator<RNASequence> redunIter = nonRedundantSequences.listIterator();
			nameRedund = new PrintWriter("Database_originale_names");
			taxaRedund = new PrintWriter("Database_originale_taxa");
			while (redunIter.hasNext()) {
				RNASequence currentSequence = redunIter.next();
				redundNames.add(currentSequence.getTaxonomy().getID());
				redundTaxa.add(currentSequence.getTaxonomy().getEMBLTaxonomyPath());
			}
			Collections.sort(redundNames);
			Collections.sort(redundTaxa);
			ListIterator<String> nameIter = redundNames.listIterator();
			while (nameIter.hasNext()) {
				nameRedund.println(nameIter.next());
			}
			nameRedund.close();
			nameIter = redundTaxa.listIterator();
			while (nameIter.hasNext()) {
				taxaRedund.println(nameIter.next());
			}
			taxaRedund.close();
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}
	}
	
/**	public void subsampleSedimentReads(String flowFile, String fastaFile, String numberSubSamples) {
		try {
			Integer numSubSampl = new Integer(numberSubSamples);
			Scanner scanFlow = new Scanner(new FileInputStream(flowFile));
			scanFlow.nextLine();
			while (scanFlow.hasNextLine()) {
				String flowString = scanFlow.nextLine();
				StringTokenizer stringToken = new StringTokenizer(flowString);
				Flow newFlow = new Flow();
				newFlow.setId(stringToken.nextToken());
				Integer cf = new Integer(stringToken.nextToken());
				newFlow.setCleanFlows(cf.intValue());
				int counter = 0;
				while (stringToken.hasMoreTokens()) {
				Double fv = new Double(stringToken.nextToken());
					newFlow.addValue(counter, fv.doubleValue());
					counter++;
				}
				subSampleFlowSequences.add(newFlow);
			}
			scanFlow.close();
			FileInputStream scanFasta = new FileInputStream(fastaFile);
			FastaReader<DNASequence, NucleotideCompound> fastaReader = 
				new FastaReader<DNASequence, NucleotideCompound>(scanFasta, new GenericFastaHeaderParser<DNASequence, NucleotideCompound>(), 
																	new DNASequenceCreator(AmbiguityDNACompoundSet.getDNACompoundSet()));
			LinkedHashMap<String, DNASequence> fastaSeq = new LinkedHashMap<String, DNASequence>();
			fastaSeq = fastaReader.process();
			subSampleFastaSequences = fastaSeq.values();
			scanFasta.close();
			int numbPerSample = subSampleFlowSequences.size() / numSubSampl;
			for (int i = 1; i < numSubSampl; i++) {
				String sampleFlowFile = "";
				String sampleFastaFile = "";
				if (i < 10) {
					sampleFlowFile = "Sample0" + i + ".dat";
					sampleFastaFile = "Sample0" + i + ".fa";
				}
				else {
					sampleFlowFile = "Sample" + i + ".dat";
					sampleFastaFile = "Sample" + i + ".fa";
				}
				PrintWriter writeSampleFlow = new PrintWriter(sampleFlowFile);
				PrintWriter writeSampleFasta = new PrintWriter(sampleFastaFile);
				writeSampleFlow.println(numbPerSample + " " + 720);
				Random randGen = new Random();
				for (int j = 1; j <= numbPerSample; j++) {
					int index = randGen.nextInt(subSampleFlowSequences.size());
					Flow currentFlow = subSampleFlowSequences.get(index);
					writeSampleFlow.print(currentFlow.getId());
					writeSampleFlow.print(" " + currentFlow.getCleanFlows());
					for (int k = 0; k < currentFlow.getSize(); k++) {
						writeSampleFlow.print(" " + currentFlow.getValue(k));
					}
					writeSampleFlow.println();
					Iterator<DNASequence> dnaIter = subSampleFastaSequences.iterator();
					boolean found = false;
					while (!found && dnaIter.hasNext()) {
						DNASequence currentSequence = dnaIter.next();
						if (currentSequence.getAccession().getID().equalsIgnoreCase(currentFlow.getId())) {
							found = true;
							writeSampleFasta.println(">" + currentSequence.getAccession().getID());
							writeSampleFasta.println(currentSequence.getSequenceAsString());
						}
					}
					if (!found) {
						System.out.println("Sequence not found in the fasta file!");
						System.exit(0);
					}
					subSampleFlowSequences.remove(index);
				}
				writeSampleFlow.close();
				writeSampleFasta.close();
			}
			String sampleFlowFile = "Sample" + numSubSampl + ".dat";
			String sampleFastaFile = "Sample" + numSubSampl + ".fa";
			PrintWriter writeSampleFlow = new PrintWriter(sampleFlowFile);
			PrintWriter writeSampleFasta = new PrintWriter(sampleFastaFile);
			writeSampleFlow.println(subSampleFlowSequences.size() + " " + 720);
			ListIterator<Flow> flowIter = subSampleFlowSequences.listIterator();
			
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		} catch (Exception e) {
			e.printStackTrace();
		}
	}
**/
//--------------------------------------------------------------------------------------------------------------------------------------------
	public static void main(String[] args) throws Exception {
		OpenFasta running = new OpenFasta();
		running.initialize(args);
	}
}